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1  Statement  of  the  Problem  Studied 


This  study  has  examined  how  to  inspect  Resident  Space  Objects  (RSOs) 
using  other  Observer  satellites,  all  of  which  are  in  fixed  orbits.  Opti¬ 
mal  methods  for  determining  which  Observers  should  characterize  each 
RSO  at  every  sample  time  have  been  developed.  The  resulting  inter¬ 
satellite  measurements  are  typically  at  a  sufficient  range  that  an  affine 
model  of  the  imaging  is  highly  accurate,  so  affine  transformations  have 
also  been  studied.  Once  these  measurements  are  taken,  they  must  then 
be  transmitted  to  Earth  through  a  communications  network,  thus  con¬ 
gestion  control  and  routing  algorithms  especially  suited  for  this  Space 
Situational  Awareness  (SSA)  task  have  been  derived.  This  section  will 
state  the  problem  studied,  with  further  details  in  the  next  section. 

This  study  developed  methods  for  regularly  characterizing  satellites 
at  high  resolution  and  low  recurring  cost  using  satellites  whose  com¬ 
bined  orbits  collectively  have  excellent  views  of  space  objects.  Humans 
and  autonomous  spacecraft  effectively  collaborate  by  allowing  humans 
to  modify  and  direct  both  overall  objectives  and  short  term  goals  while 
using  autonomous  planning  techniques  to  manage  complex,  yet  well- 
defined  and  predictable  orbital  events  and  sensor  allocations.  This 
allows  humans  to  provide  input  where  they  are  strong,  such  as  un¬ 
derstanding  high  level  goals  and  reasoning  amidst  uncertainty.  Simul¬ 
taneously,  it  allows  autonomous  optimization  and  planning  techniques 
to  operate  in  domains  where  they  are  strong,  in  this  case  when  many 
spacecraft  must  be  jointly  controlled  over  space  and  time  to  achieve  a 
clearly  defined  overall  final  result.  This  study  formulates  observation  of 
multiple  satellites  by  a  set  of  Observer  satellites  first  as  a  binary  integer 
programming  problem.  Next,  methods  are  found  to  relax  the  formu¬ 
lation  to  become  a  convex  linear  program.  Thus  global  optimality, 
high  speed,  and  convergence  to  the  optimum  are  assured.  Real-time 
human  input  is  allowed  by  modification  of  the  optimization  weights. 
Simulation  results  confirm  the  viability  of  the  technique  on  a  variety 
of  orbits,  including  both  low  Earth  and  geosynchronous  cases.  In  ad¬ 
dition,  the  simulations  indicate  that  the  binary  approximation  of  the 
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convex  relaxation  achieves  nearly  optimal  performance.  The  method  is 
fully  general,  encompassing  3-D  and  elliptical  orbits. 

Once  measurements  are  taken,  they  must  be  transmitted  to  Earth, 
so  communications  networks  especially  suited  for  this  SSA  task  have 
been  developed.  A  stable  congestion  controller  using  data-driven,  safe 
switching  control  theory  is  developed  that  improves  the  dynamic  per¬ 
formance  of  the  satellite  networks.  Robustness  is  improved  by  allowing 
the  switching  among  members  of  a  family  of  candidate  control  laws, 
rather  than  on  using  a  fixed  controller  based  on  a  priori  assumptions 
about  the  network  dynamics.  The  stable  region  of  the  Proportional- 
Integral  (PI)  parameters  for  a  nominal  model  is  explored  first.  Then, 
a  PI  controller,  whose  parameters  are  adaptively  tuned  by  switching 
among  members  of  a  given  candidate  set  using  observed  plant  data,  is 
presented  and  compared  with  several  standard  AQM  policy  examples. 
A  new  cost  detectable  switching  law  with  the  interval  cost  function 
switching  algorithm  is  developed,  which  improves  the  performance  and 
saves  the  computational  cost;  it  is  then  compared  with  a  law  com¬ 
monly  used  in  the  switching  control  literature.  Finite-gain  stability  of 
the  system  is  proven. 

In  addition  to  congestion  control,  an  efficient  routing  algorithm  plays 
a  key  role  in  optimizing  network  resources  in  a  LEO  satellite  net¬ 
work.  A  new  routing  system  based  on  the  Cross  Entropy  optimiza¬ 
tion  method  has  been  developed  in  this  study.  The  novel  on-demand 
routing  method,  named  Cross  Entropy  Accelerated  Ant  Routing  Sys¬ 
tem  (CEAARS)  is  applied  in  simulations  to  the  regular  constellation 
LEO  satellite  networks.  Simulations  on  an  Iridium-like  satellite  net¬ 
work  compare  the  proposed  CEAARS  algorithm  with  two  standard 
approaches  to  adaptive  routing  protocols  on  the  Internet:  Distance- 
Vector  (DV)  and  Link-State  (LS),  as  well  as  with  the  original  Cross 
Entropy  Ant  Routing  System  (CEARS).  DV  algorithms  are  based  on 
the  distributed  Bellman  Ford  algorithm,  and  LS  algorithms  are  an  im¬ 
plementation  of  Dijkstra’s  single  source  shortest  path  algorithm.  The 
results  show  that  CEAARS  not  only  remarkably  improves  the  conver¬ 
gence  speed  of  achieving  optimal  or  suboptimal  paths,  but  also  reduces 
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the  number  of  overhead  ants  (management  packets). 

From  a  single  constellation  of  satellites,  several  layers  are  needed 
to  transmit  the  information  to  Earth.  Thus  a  multi-layered  satellite 
network  has  also  been  developed.  Due  to  the  dynamic  topology  of  the 
network,  large  propagation  delays  and  bulk  data  transfers,  significant 
delays  and  loss  of  data  can  be  encountered  with  traditional  TCP  con¬ 
gestion  control  schemes.  To  overcome  these  drawbacks,  network  snap¬ 
shots  consisting  of  layered  dynamic  clusters  of  satellites  are  formed 
with  modified  TCP  congestion  control  schemes  such  that  the  delay  and 
packet  losses  are  at  a  minimum.  Using  our  new  convex  optimization 
formulation,  optimal  incoming  data  rates  are  computed  at  each  satel¬ 
lite  irrespective  of  their  layers.  They  are  computed  in  each  snapshot  to 
keep  the  delay  and  packet  loss  at  a  minimum.  During  the  process  of 
determining  optimal  incoming  data  rates  a  collection  of  satellites  are 
also  identified  to  perform  a  particular  space  situational  awareness  task. 

The  resulting  inter-satellite  measurements  are  typically  at  a  suffi¬ 
cient  range  that  an  affine  model  of  the  imaging  is  highly  accurate, 
so  affine  transformations  and  related  approximation  methods  have  also 
been  studied.  The  celebrated  Paley-Wiener  theorem  naturally  identifies 
the  spaces  of  bandlimited  functions  with  subspaces  of  entire  functions 
of  exponential  type.  This  study  shows  that  these  spaces  remain  invari¬ 
ant  only  under  composition  with  affine  maps.  After  some  motivation 
demonstrating  the  importance  of  characterization  of  range  spaces  aris¬ 
ing  from  the  action  of  more  general  composition  operators  on  the  spaces 
of  bandlimited  functions,  the  subspaces  of  L2(R)  generated  by  these  ac¬ 
tions  are  identified.  Extension  of  these  theorems  where  Paley-Wiener 
spaces  are  replaced  by  the  deBranges-Rovnyak  spaces  are  given. 


2  Summary  of  the  Most  Important  Results 

A  more  detailed  summary  of  this  study  will  now  be  presented  in  the 
following  subsection. 
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2.1  Optimal  Observation  of  Satellites  Using  Combined  Mea¬ 
surements  from  Many  Networked  Observers  Coopera¬ 
tively  Controlled  by  Human  and  Autonomous  Means 


New  methods  for  regularly  characterizing  satellites  at  high  resolution 
and  low  recurring  cost  using  satellites  whose  combined  orbits  collec¬ 
tively  have  excellent  views  of  space  objects  have  been  developed.  Hu¬ 
mans  and  autonomous  spacecraft  effectively  collaborate  by  allowing 
humans  to  modify  and  direct  both  overall  objectives  and  short  term 
goals  while  using  autonomous  planning  techniques  to  manage  complex, 
yet  well-defined  and  predictable  orbital  events  and  sensor  allocations. 
This  allows  humans  to  provide  input  where  they  are  strong,  such  as  un¬ 
derstanding  high  level  goals  and  reasoning  amidst  uncertainty.  Simul¬ 
taneously,  it  allows  autonomous  optimization  and  planning  techniques 
to  operate  in  domains  where  they  are  strong,  in  this  case  when  many 
spacecraft  must  be  jointly  controlled  over  space  and  time  to  achieve 
a  clearly  defined  overall  final  result.  Characterization  of  satellites  in 
geosynchronous  orbit  is  especially  important,  as  they  are  too  far  away 
for  high  resolution  observation  from  Earth.  This  report  formulates  ob¬ 
servation  of  multiple  satellites  by  a  set  of  Observer  satellites  first  as 
a  binary  integer  programming  problem.  Next,  methods  are  found  to 
relax  the  formulation  to  become  a  convex  linear  program.  Thus  global 
optimality,  high  speed,  and  convergence  to  the  optimum  are  assured. 
Real-time  human  input  is  allowed  by  modification  of  the  optimization 
weights.  For  instance,  if  a  human  decides  that  a  particular  object  re¬ 
quires  extra  observation  immediately,  then  this  paper  shows  how  this 
short  term  goal  can  be  achieved  by  weight  modification.  Similarly,  if 
a  human  operator  desires  that  an  object  be  observed  more  persistently 
and  accurately  than  other  objects,  then  this  paper  develops  weighting 
methods  to  achieve  this  long  term  goal.  Simulation  results  confirm  the 
viability  of  the  technique  on  a  variety  of  orbits,  including  both  low 
Earth  and  geosynchronous  cases.  In  addition,  the  simulations  indicate 
that  the  binary  approximation  of  the  convex  relaxation  achieves  nearly 
optimal  performance.  The  method  is  fully  general,  encompassing  3-D 
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and  elliptical  orbits. 

In  order  to  automatically  plan  these  satellite  to  satellite  observa¬ 
tions,  we  have  developed  a  new  relaxed  Integer  Linear  Programming 
(relaxed-ILP)  algorithm  which  can  rapidly  generate  optimal  plans.  Fig. 
1-4  depict  a  larger  scale  simulation  (n=40  RSOs,  m=10  Observers, 
N=201  samples)  of  nearly  Geostationary  orbits.  The  number  of  RSOs 
and  Observers  is  representative  of  a  true  GEO  characterization  task, 
as  there  are  approximately  400  active  GEO  satellites  currently,  thus 
this  scenario  would  characterize  the  most  important  10  %  of  the  GEO 
satellites.  The  total  simulation  time  covers  two  orbits.  The  RSOs  are 
rotating  the  same  direction  as  the  Earth  (as  is  customary),  while  the 
Observers  rotate  in  the  opposite  direction  so  they  will  pass  nearby  all  of 
the  RSOs  twice  per  orbit.  Fig.  1  illustrates  the  first  sample  time-due 
to  the  high  GEO  altitude,  occlusion  by  the  Earth  is  not  an  issue,  and 
all  RSOs  are  visible  and  are  actively  being  sensed  by  some  Observer. 
Each  Observer  spends  most  of  the  sample  time  on  a  single  RSO  but 
also  senses  other  RSOs  for  a  small  portion  of  the  sample  time.  Note 
that  this  result  is  impractical,  as  the  Observers  are  attempting  to  slew 
towards  too  many  RSOs  per  sample,  and  many  RSOs  are  far  away. 
These  problems  can  be  partially  corrected  through  a  binary  approxi¬ 
mation  which  forces  each  Observer  to  view  only  one  RSO  per  sample, 
but  slewing  issues  become  more  pronounced  as  the  slewing  rates  be¬ 
come  high,  as  in  the  practical  retro-GEO  case.  Fig.  2  depicts  the 
results  halfway  through  the  simulation  (after  a  single  orbit  has  been 
completed).  At  this  point,  there  is  still  some  marked  variability  in  the 
quality  of  each  RSO’s  characterization  (as  indicated  by  the  size  of  the  X 
marks).  Fig.  3  illustrates  the  final  results-many  observations  are  tak¬ 
ing  place  simultaneously,  and  the  final  results  are  excellent-all  RSOs 
have  been  observed  nearly  equally. 

Fig.  4  illustrates  the  results  of  the  binary  approximation,  which 
concentrates  all  sensing  on  a  single  RSO  per  Observer  during  each 
sample  time  to  crudely  reduce  slewing.  In  this  simulation,  it  displays 
high  performance  that  is  very  close  to  optimal  because  all  RSOs  are  of 
very  nearly  equal  size,  which  indicates  that  they  have  been  observed 


at  nearly  equal  quality.  Note  that  during  the  last  sample,  each  Ob¬ 
server  sensed  either  the  first  or  second  closest  RSO,  depending  upon 
the  optimal  plan. 

To  improve  the  binary  approximation,  the  relaxed-ILP  algorithm 
just  illustrated  has  been  combined  with  a  Greedy  optimization  which 
finds  the  optimal  binary  threshold.  This  new  ILP-Greedy  algorithm  is 
a  amalgam  of  the  relaxed  ILP  and  a  Batch  Greedy  algorithm,  and  it  is 
based  on  the  block  diagram  illustrated  in  Figure  5.  As  shown,  first  the 
ILP  solution  is  found  using  the  ILP  relaxation  method,  then  a  binary 
approximation  will  be  done  by  an  optional  threshold  t  >  0.5. 

After  the  binary  approximation  process,  the  plan  will  not  have  used 
all  Observers  at  all  sample  times  because  a  high  threshold  sets  all  ele¬ 
ments  in  those  columns  to  zero.  In  the  last  step  of  this  mixed-algorithm, 
the  batch  greedy  algorithm  will  be  employed  to  fill  these  zero  columns. 
Note  that  when  the  threshold,  t,  equals  zeros  then  a  high  performance 
but  impractical  pure  linear  program  is  used,  while  when  t  =  1,  the 
ILP-Greedy  algorithm  is  nothing  but  a  batch  greedy  algorithm.  As  the 
Batch  Greedy  algorithm  is  fast,  one  can  run  the  algorithm  with  different 
values  of  t  and  choose  the  best  planning  among  those  to  get  a  practical, 
but  slightly  sub-optimal  plan.  Simulation  results  show  that  the  mixed 
algorithm  with  optimal  threshold  values  always  performs  better  than 
either  the  relaxed-ILP  or  the  Batch  Greedy  algorithms. 

Figures  6  and  7  illustrate  typical  results  of  all  discussed  planning 
methods  for  two  different  planning  problems.  All  satellites  in  the  prob¬ 
lem  related  to  Figure  6  are  arbitrarily  chosen  within  the  Low  Earth  Or¬ 
bit  (LEO),  but  all  satellites  related  to  the  other  figure  orbit  in  Geosyn¬ 
chronous  orbit  (GEO). 

The  planning  problem  in  LEO  is  a  comparatively  small  problem 
(n  =  5  RSOs,  m  =  3  Observers,  N  =  500  Samples  and  L  =  2  side  of 
interests  — >  5x3x500  =  7500  unknowns),  while  the  application  in  GEO 
is  a  much  bigger  planning  problem  (n  =  40  RSOs,m  =  8  Observers, 
N  =  901  Samples  and  L  =  2  side  of  interests  -4  40  x  8  x  901  = 
288320  unknowns).  In  both  Figures,  black  dotted,  green  dashed  and  red 
dashed-dot  lines  imply  the  upper  bound  obtained  from  ILP-Relaxation 
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Orbit  of  Satellites  Around  Earth  for  sample  1/201 


x  104 


Figure  1:  GEO  Orbital  Simulation  for  the  first  sample  time.  Observers  are  indicated 
by  green  circles  and  RSOs  by  red  x’s.  The  width  of  a  line  connecting  an  Observer  to  an 
RSO  is  proportional  to  the  percentage  of  that  sample  time  during  which  the  Observer 
collects  data  from  the  corresponding  RSO,  and  the  size  of  the  RSO  is  proportional 
to  the  total  quality  of  its  observation  at  the  current  sample  time 
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Orbit  of  Satellites  Around  Earth  for  sample  100/201 
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Figure  2:  GEO  Orbital  Simulation  after  one  orbit. 
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Orbit  of  Satellites  Around  Earth  for  sample  201/201 
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Figure  3:  GEO  Orbital  simulation  for  the  last  sample  time. 
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Orbit  of  Satellites  Around  Earth  for  sample  201/201 


x  104 


x  104 


Figure  4:  GEO  Binary  orbital  simulation  depicted  at  the  last  sample  time.  At  this 
stage  the  result  becomes  practical,  as  all  Observers  are  measuring  nearby  RSOs,  with 
only  one  RSO  viewed  each  sample  time  by  each  Observer. 
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Figure  5:  The  block  diagram  of  ILP-Greedy  method. 


Performance  Comparision  :  All  Planning  Algorithms 


Figure  6:  Performance  comparison  between  different  planning  algorithms  for  LEO 
satellites. 


method,  solution  of  Batch  Greedy  algorithm  and  solution  of  the  Binary 
approximation  of  relaxed  ILP,  respectively.  The  x  axis  is  used  to  show 
different  values  of  the  binary  approximation  threshold  and  the  y  axis 
points  out  the  minimum  observation  quality  of  all  RSOs  in  all  sides 
of  interests.  A  red  circle  marker  is  used  to  show  the  best  solution  of 
ILP-Greedy  method  for  the  different  threshold  values. 
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Performance  Comparision  :  All  Planning  Algorithms 


Figure  7:  Performance  comparison  between  different  planning  algorithms  for  GEO 
satellites. 


The  most  noticeable  feature  of  the  results  is  that  the  solution  of  ILP- 
Greedy  algorithm  is  closer  to  the  upper  bound  than  the  other  methods, 
and  it  simply  means  the  ILP-Greedy  offers  a  better  solution.  Moreover, 
in  the  LEO  problem,  the  ILP  method  performs  better  than  the  batch 
Greedy  algorithm,  but  in  the  GEO  problem  the  converse  applies,  and 
it  confirms  that  neither  of  these  two  approximate  methods  is  preferable 
in  general.  A  very  powerful  feature  of  this  overall  approach  is  that 
an  upper  bound  on  the  gap  between  the  found  sub-optimal  solution 
and  the  unknown  optimal  solution  is  available.  The  ILP-relaxation 
algorithm  provides  an  optimal  but  physically  unrealizable  solution.  So 
if  performance  approaches  it  (  as  it  does  in  both  Figures  6  and  7),  then 
the  sub-optimal  solution  is  very  nearly  optimal. 

Figure  8  depicts  a  snapshot  of  an  observation  planning  for  the  LEO 
problem  at  sample  time  k  =  158.  Observers  are  indicated  by  numbered 
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Figure  8:  A  snapshot  ILP-Greedy  planning  in  LEO  orbit  at  sample  time  K  =  158. 

red  circles  and  RSOs  by  numbered  blue  squares.  All  satellites  are  mov¬ 
ing  in  six  different  orbits,  3  RSO  orbits  and  3  Observer  orbits.  Black 
lines  are  utilized  to  indicate  that  a  RSO  is  inspected  by  the  connecting 
observer  at  this  time  instance.  Qsum  the  denotes  cumulative  observa¬ 
tion  quality  matrix,  and  it  shows  the  observation  qualities  of  all  RSOs 
in  every  individual  side  of  interest.  Cartesian  coordinate  axes  x,  y  and 
z  are  shown  in  red  ,  blue  and  light  green  respectively.  The  direction  of 
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Rso3 


Observed  Observed  Observer3  Observer4 

Figure  9:  Relative  orbits  of  Observers  with  respect  to  RSO 


sun  light  is  from  positive  to  negative  side  of  the  x  axis. 

All  possible  RSOs  images  that  can  be  taken  from  all  Observers  at 
the  k  =  158  time  instance  are  gathered  in  Figure  8.  Every  observer  has 
five  options  (RSOs)  to  choose  for  improving  the  observation  quality. 
The  following  items  explain  the  reasoning  of  choosing  a  RSO  for  every 
individual  observer.  Note  that  this  reasoning  is  built  into  the  auto¬ 
matic  planning  system  (by  minimizing  performance  indices,  not  logical 
inference). 
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Observer  1 


According  to  Figure  8  RSO  3  and  5  are  occluded  by  the  Earth  so  the 
observation  quality  is  zero.  Although  the  distance  between  this  observer 
and  RSO  1  is  the  shortest,  due  to  sun  light  direction  the  observer  only 
can  inspect  the  dark  side  of  it.  Even  though  the  value  of  observation 
quality  for  RSO  4  is  greater  than  that  of  RSO  2,  the  ILP-Greedy  scheme 
chooses  RSO  2  because  it  knows  other  good  opportunities  for  observing 
RSO  4  will  occur,  and  the  long  term  planning  is  possible.  In  contrast, 
the  Batch-Greedy  process  lacks  this  long  term  perspective,  so  the  Batch 
Greedy  Algorithm  would  choose  RSO  1. 

Observer  2 

RSO  1  ,  2  and  4  are  occluded  by  the  Earth.  Despite  the  fact  that  RSO 
3  and  5  are  relatively  close  to  observer  1,  the  deep  darkness  of  that  area 
makes  those  images  worthless. 

Observer  3 

Again  Earth  occlusion  causes  Observer  3  not  to  see  RSO  3  and  5,  and 
RSO  2  and  4  Images  are  blurred  because  of  the  large  distance.  So,  the 
ILP-Greedy  algorithm  chooses  RSO  1  for  inspection.  In  this  case,  the 
Greedy  algorithm  also  chooses  RSO  1. 

Relative  Orbit  Analysis 

Figure  9  depicts  relative  orbits  of  all  the  Observers  with  respect  to  the 
RSOs  in  another  application  within  LEO,  and  it  shows  that  the  position 
of  each  Observer  with  respect  to  each  RSO  satellite,  in  a  frame  attached 
to,  and  rotating  with,  the  RSO  satellite  (Hill’s  frame).  Black  lines 
imply  all  observations  made  by  observers  in  different  sample  times.  By 
analyzing  relative  orbits  of  each  RSO  along  with  planning,  it  is  easier 
to  evaluate  the  performance  of  the  proposed  method. 
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2.2  Stable  switching  congestion  control  for  satellite  TCP / AQM 
networks 

Once  measurements  have  been  taken,  they  must  be  eventually  trans¬ 
mitted  to  Earth.  Methods  for  doing  so  that  are  suited  for  this  SSA 
scenario  are  now  outlined. 

A  stable  congestion  controller  using  data-driven,  safe  switching  con¬ 
trol  theory  is  developed  that  improves  the  dynamic  performance  of  the 
satellite  TCP/ AQM  networks  [1],  [2].  Robustness  is  improved  by  allow¬ 
ing  the  switching  among  members  of  a  family  of  candidate  control  laws, 
rather  than  on  using  a  fixed  controller  based  on  a  priori  assumptions 
about  the  network  dynamics.  The  stable  region  of  the  Proportional- 
Integral  (PI)  parameters  for  a  nominal  model  is  explored  first.  Then, 
a  PI  controller,  whose  parameters  are  adaptively  tuned  by  switching 
among  members  of  a  given  candidate  set  using  observed  plant  data,  is 
presented  and  compared  with  several  standard  AQM  policy  examples, 
such  as  Random  Early  Detection  (RED)  and  fixed  PI  control.  A  new 
cost  detectable  switching  law  with  the  interval  cost  function  switching 
algorithm  is  developed,  which  improves  the  performance  and  saves  the 
computational  cost;  it  is  then  compared  with  a  law  commonly  used  in 
the  switching  control  literature.  Finite-gain  stability  of  the  system  is 
proven.  More  details  follow  in  the  text  below. 

Congestion  on  the  satellite  networks  is  one  of  the  major  communica¬ 
tion  problems.  Transmission  Control  Protocol  (TCP)  congestion  con¬ 
trol  algorithms  use  packet-loss  and  packet-delay  measurements  respec¬ 
tively  to  detect  congestion.  The  current  standard  in  controlling  end-to- 
end  congestion  in  the  Internet  is  the  so-called  Active  Queue  Manage¬ 
ment  (AQM);  its  functionality  is  based  on  sensing  impending  congestion 
before  it  occurs  and  providing  feedback  information  to  senders  by  either 
dropping  or  marking  packets,  so  that  congestion,  causing  a  significant 
degradation  in  network  performance,  can  be  avoided. 

Prior  developments  in  the  robust  congestion  control  considered  only 
a  single  specific  model  of  the  TCP/AQM  dynamics,  which  virtually 
reduced  the  effectiveness  of  the  control  scheme  for  all  other  than  the 
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nominal  operating  point  conditions.  The  final  results  of  this  work  build 
on  our  previous  results  in  the  Safe  Switching  Control  (SSC)  of  uncer¬ 
tain  TCP  systems  [2],  and  develop  a  new  cost  detectable  switching  law 
with  the  interval  cost  function  switching  algorithm,  which  improves  the 
performance  while  saving  the  computational  cost;  comparisons  with  a 
law  commonly  used  in  the  switching  control  literature  are  provided. 
Finite-gain  stability  of  the  system  is  proven.  A  fuzzy  logic  PI  con¬ 
troller  is  proposed  as  a  special  candidate  to  achieve  good  performance 
at  all  nominal  points  with  the  available  set  of  candidate  controllers 
(details  in  [3]  and  [4]).  The  new  AQM  algorithm  is  applicable  to  long 
fat  networks  and  in  particular  multi-layer  satellite  networks.  Through 
Matlab  and  ns- 2  simulations,  we  compare  the  performance  of  the  new 
SSC  scheme  with  that  of  Random  Early  Detection  (RED),  adaptive 
RED,  and  fixed  PI  controller.  From  the  simulation  results,  it  is  shown 
that  the  proposed  SSC  scheme  can  adaptively  deal  with  the  change  of 
the  number  of  connections  and  heterogeneous  delays. 

To  design  congestion  controllers,  large  scale  networks  are  often  sim¬ 
plified  according  to  a  ‘dumbbell’  topology,  with  the  network  consisting 
of  n  senders,  one  bottleneck  router  and  one  receiver,  which  in  a  clus¬ 
ter  form  of  satellite  networks  correspond  to  normal  satellite  nodes,  a 
cluster-head  of  one  cluster  and  a  cluster-head  of  another  (adjacent) 
cluster,  respectively.  In  this  topology,  we  assume  that  N  TCP  connec¬ 
tions  represent  homogeneous  and  long  lived  flows.  A  fluid-flow  model 
describing  the  behavior  of  the  average  values  of  key  network  variables 
(window  size  in  senders,  Round-Trip  Time  (RTT),  and  the  queue  size 
in  the  bottleneck  router)  is  given  by  the  following  coupled,  nonlinear 
differential  equations, 


W(t) 


1 

W) 


W(t)W(t  -  R(t )) 
2 R(t  -  R(t )) 


■p{t  -  R(t )) 


<i(t)  = 


m 


N(t )  -  C 


(1) 

(2) 


where  R{t)  =  ^=r  +  Tp  is  the  RTT,  IV  is  the  average  TCP  window  size,  p 
is  the  dropping  probability  of  a  packet,  N  is  the  number  of  connections 
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or  TCP  sessions,  C  is  the  transmission  capacity  of  the  router,  q  is  the 
average  queue  length  of  the  router  buffer,  and  Tp  is  the  propagation 
delay.  For  the  design  of  a  nominal  plant  model,  the  system  is  linearized 
around  the  equilibrium  point  (Wq,  qo,po). 

In  the  prior  works,  some  improvements  in  the  performance  (queue 
utilization,  packet  drop)  have  been  achieved  using  adaptive  RED  and 
PI  schemes;  however,  all  these  methods  are  based  on  the  linearization 
model  at  one  equilibrium  point  only,  and  thus  are  not  effective  in  the 
multi-layer  satellite  networks,  where  both  N  and  R  are  dynamically 
varying. 

Here  we  exploit  the  results  on  time-delay  systems  to  characterize 
the  stable  region  of  kp  and  kj  gains.  To  improve  the  performance  in 
the  dynamically  changing  operating  conditions,  we  consider  first  the 
PI  stabilization  of  first  order  systems  with  time  delay  and  extend  to 
the  case  of  TCP  AQM.  The  characteristic  equation  of  dynamic  models 
with  time  delays  has  the  form: 

6(s)  =  d(s)  +  e~LlSni(s)  +  e~L2Sn2(s)  H - P  e~LmSnm(s)  (3) 

where  d(s)  and  rq(s)  for  i  =  1,2,  •••  ,  m  are  polynomials  with  real 
coefficients;  Lt  is  the  time  delay;  deg[d(s)]  =  q  and  deg[rii(s)]  <  q; 
0  <  L\  <  L2  <  ■  ■  ■  <  Lm.  Instead  of  (3)  we  can  consider  the  quasi¬ 
polynomial  3*  as 

6*{s)  =  eLmS5(s ) 

=  eLmSd(s )  +  e^Lm~L^sni(s)  + 

e(Lm-L2)sn2(g)  - p  nm(s )  (4) 

To  guarantee  the  stability  of  the  system,  all  zeros  of  the  characteristic 
equation  (3)  (and  thus,  all  zeros  of  (4))  must  be  in  the  open  left  half 
plane.  The  following  theorem  gives  a  necessary  and  sufficient  condition 
for  the  stability  of  (4). 

Theorem  1.  Rewrite  3*(s)  in  (4)  as 

s*(ju)  =  Sr(u)  +  jSi(u)  (5) 
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where  5r( u)  and  Sfu)  respectively  represent  the  real  and  imaginary  parts 
ofS*(jw).  Then  5*(s)  is  stable  if  and  only  if: 

1:  5r( u)  and  dfuj)  have  only  simple,  real  roots  and  these  interlace. 

2:  ^'(cJo)^r(^o)  —  >  0  for  some  uo  in  (— oo,  oo) 

where  S'r(cuo)  and  8'(cu o)  denote  the  first  derivative  with  respect  to  u  of 

5r(cu)  and  Sfu),  respectively. 

The  following  theorem  verifies  that  the  5r(u)  and  Sfcv)  have  only 
real  roots. 


Theorem  2.  Let  M  and  N  denote  the  highest  powers  of  s  and  es , 
respectively  in  8*(s).  Let  rj  be  an  appropriate  constant  such  that  the 
coefficients  of  terms  of  highest  degree  in  5r((j)  and  Sfui)  do  not  vanish 
at  uj  =  rj.  Then  for  the  equations  5r{u)  =  0  or  8j(u)  =  0  to  have  only 
real  roots,  it  is  necessary  and  sufficient  that  in  each  of  the  intervals 
—2ln  +  rj  <  cu  <  21tt  +  rj  with  l  =  Iq,Iq  +  1,  lo  +  2  •  •  • ,  5r(uj)  and  5i(u) 
have  exactly  LIN  +  M  real  roots  for  a  sufficiently  large  Iq. 


PI  Control  for  the  TCP  Plant  Dynamics 

The  linearized  TCP  plant  can  be  rewritten  as  the  following  transfer 
function: 


P(s)  = 


K 


— Ls 


(6) 


s2  +  a\s  +  ao 

where  K  is  the  static  gain  of  the  plant,  L  is  the  time  delay,  and  K ,  L ,  ao 
and  a\  are  all  nonnegative  parameters.  Our  objective  is  to  analytically 
determine  the  set  of  all  stabilizing  (kp,  ki )  values  for  the  given  plant. 
The  characteristic  equation  of  the  closed  loop  system  is: 


S(s)  —  I\  {kj  T  kps)e  ^ s  T  (s^  T  oqs  T  ag)s  (7) 


To  use  Theorem  1  to  hnd  the  set  of  stabilizing  PI  controllers,  we  deduce 
the  quasi-polynomial  8*(s),  he., 

8*(s)  =  5(s)eLs  =  K{ki  +  kps )  +  (s2  +  a\s  +  ao)seis  (8) 


by  replacing  s  by  ju,  we  get, 

8*(ju)  =  8r{uj)  +j8i(u) 


(9) 
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By  sweeping  over  all  real  kp  and  solving  a  constant  gain  stabilization 
problem  at  each  stage,  we  can  determine  the  region  in  the  (kp,  kf) 
parameter  space  for  which  the  closed-loop  system  is  stable.  We  will 
assume  kp  >  0,  kj  >  0.  Then,  the  range  of  kp  over  which  the  sweeping 
needs  to  be  carried  out  can  be  narrowed  down  by  using  the  following 
theorem: 


Theorem  3.  Under  the  above  assumptions  on  K,  L,  ao  and  a\,  the 
range  of  kp  values  for  which  a  solution  exists  to  the  PI  stabilization 
problem  of  a  given  open-loop  stable  plant  with  transfer  function  P(s )  as 
in  equation  (6)  is  given  by 

1  ca  of 

0  <  kP  <  —  (ai- sin(a)  -  cos(o)(a0  -  t^))  (10) 

K  L  lv 

Where  a  is  the  solution  of  the  equation 


tan(a) 


in  the  interval  [0, 7r  . 


cx( 2  T  a\L) 
a1  —  ao L2  ~  O'lL 


ai) 


Proof.  By  substituting  z  =  Luj ,  the  real  and  imaginary  parts  of  S*(s) 
can  be  rewritten  as, 

z 

—  sin 


(*))  (12) 


z  z z 

Si(z)  =  — ( Kkp  +  (a0  -  J?)  cos (z)  -  ai 


z^  z  z^ 

Sr(z )  =  Kki  +  (—  -  a0— )  sin(^)  -  aij^  cos  (z) 
=  K[kj  —  a(z 


(13) 


where 

9 

z  z  z 

a(z)  =  -^(sin (z)(a0  -  j^)  +  ax-  cos(2:))  (14) 

According  to  Theorem  1,  we  check  two  conditions  to  ensure  the 
stability  of  the  quasi-polynomial  S*(s). 

We  start  by  checking  condition  2  of  Theorem  1: 


E(zq)  =  S'i(zo)Sr(z0)  -  5i(zo)5'r(zo )  >  0  (15) 
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(16) 


KkP 

~T~ 


(a0  + 

,oo 

[T 


2a\Z 

~U~ 

3  z2 
'  ZA 


— ^)  sin(,z)  + 


z“ 


aijp)  cos(z) 


By  substituting  zo  =  0,  we  get  8i(z)  =  0,  8[{z)  =  —  kp+a°  and  5r(z)  = 
Kkj.  Hence, 


E(zq) 


5'i{z0)5r(zo) 


Kkp  +  ao 

Z 


Kki  >  0 


(17) 


Since  K  >  0  and  kj  >  0,  we  have  kp  >  —jk.  Together  with  previous 
assumption,  we  have  kp  >  0  and  ki  >  0. 

We  now  check  condition  1  of  Theorem  1:  the  interlacing  of  the  roots  of 
Sr(z)  and  Si(z).  From  equation  (12),  we  can  calculate  the  roots  of  the 
imaginary  part,  i.e. ,  5i(z)  =  0.  This  gives  us  the  following  equation 

2 

Z z  z 

z  =  0  or  Kkp  +  (ao  — rz)  cos(^)  —  a\—  sin(^)  =  0  (18) 

Lz  L 


From  this  we  see  that  one  root  of  the  imaginary  part  is  zq  =  0.  The 
other  roots  are  difficult  to  find  since  we  need  to  solve  (18)  analytically. 
However,  we  can  plot  the  terms  involved  in  (18)  and  graphically  ex¬ 
amine  the  nature  of  the  solution.  There  are  three  different  cases  to 
consider.  In  each  case,  the  positive  real  roots  of  (18)  will  be  denoted 

by  Zj ,  j  =  1,  2, ...,  arranged  in  increasing  order  of  magnitude. 

2 

Let  f(z)  =  Kkp  +  (ao  —  jp)  cos(z)  and  g(z)  =  ciij-sin^z) 

Case  1  :  kp  =  ku.  Where  ku  is  the  maximal  value  of  kp  so  that  the 
graphs  of  f(z)  and  g(z)  are  tangent  in  the  interval  [0, 7r] .  According  to 
the  definition  of  ku,  if  kp  =  ku,  the  curves  of  f(z)  and  g(z)  are  tangent 
at  the  point  a,  i.e. 


f  f(z  =  a)  =  g(z  =  a) 
\  f(z  =  a)  =  g'(z  =  a) 


From  equation  group  (19),  we  obtain: 

ck(2  T  cl\L ) 


tan  (ex)  = 


a2  —  clqL2  —  clij 


(19) 


(20) 
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Figure  10:  Plots  of  f(z)  and  g(z)  in  case  2 


Then,  the  parameter  ku  can  be  expressed  by  (21), 

1  a  or1 

ku  =  —  (aq- sin(cx)  -  cos(a)(a0  -  —,)) 
K  L 


(21) 


In  this  case,  we  can  graph  f(z )  and  g(z).  For  example,  when  R  = 
0.03,  N  =  80,  C  =  3750,  the  plot  is  shown  in  Fig.  10. 

Case  2:0  <  kp  <  ku.  In  this  case,  the  graph  of  f(z )  and  g(z)  is 
obtained  by  moving  down  the  line  of  f(z)  in  Fig.  10. 

Case  3  :  kp  >  ku.  In  this  case,  the  graph  of  f(z)  and  g(z)  is  obtained 
by  moving  up  the  line  of  f(z )  in  Fig.  10. 

Let  us  now  use  the  results  presented  in  Theorem  2  to  check  if  5i(z) 
has  only  simple  roots.  By  substituting  si  =  Ls  in  the  expression  for 
3*(s)  in  (8),  we  obtain: 


■  si 


■  si 


si. 


Si 


S*M  =  eSl((-^  +  ai(-f)2  +  tO  +  K{kp^  +  kj)  (22) 


L 


L 


L 


L 


We  see  that  for  the  new  quasi-polynomial  in  si,  M  =  3  and  N  =  1, 
where  M  and  N  designate  the  highest  powers  of  s  and  es  which  appear 
in  3*(s).  Next  we  choose  rj  =  |  to  satisfy  the  requirements  that  Sr(g)  ^ 
0  and  5i(rj)  ^  0.  In  case  2,  we  see  that,  for  0  <  kp  <  ku ,  Si(z)  has  four 
roots  in  the  interval  [0,  ^],  including  a  root  at  the  origin.  Since  Si(z) 
is  an  odd  function  of  z,  it  follows  that  in  the  interval  [—  ^f\  ,  5i(z) 

has  seven  roots.  Hence,  Si(z )  has  4N  +  M  =  7  roots  in  the  interval 

[~  T’  t],  01  even  [— T’  t]-  ^so  has  ^w0  rea^  r0°fs  in  each  °f  the 
interval  [2kir  +  |,  2 {k  +  1)7 r  +  |]  and  [—2 (k  +  1)7 r  +  |,  —  2kn  +  |]  for 
k  =  1,  2, ....  Thus,  for  0  <  kp  <  ku ,  5i(z)  has  4 kN  +  M  real  roots  in  the 
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interval  [— 2&7T  + |,  2/c7t  +  |].  Hence,  from  Theorem  2,  we  conclude  that 
Si(z)  has  only  real  roots  for  every  kp  in  [0,  ku\.  For  the  case  kp  >  ku,  the 
roots  of  5i(z )  are  not  real,  thereby  ruling  out  closed-loop  stability.  □ 

We  now  evaluate  Sr(z)  at  the  roots  of  the  imaginary  part  St(z).  For 
zo  =  0,  we  obtain 

Sr(z0)  =  K[kj  -  a(0)] 

=  Kkj  (23) 

For  Zj  7^  Zq  where  j  =  1,2,3,...,  using  (13)  we  obtain 

Sr(zj)  =  K[kj  —  a(zj )]  (24) 

Interlacing  the  roots  of  5r(z)  and  5i(z)  is  equivalent  to  Sr(zo)  >  0  (since 
ki  >  0),  Sr(zi)  <  0,  Sr(z2)  >  0,  5r (^3)  <  0,  and  so  on.  Using  this  fact 
and  (23,  24)  we  obtain 

5r(zi)  <0  =^>  kj  <  a(zi) 

Sr(z 2)  >0  =4>  kj  >  a(z2 ) 
dr(z3)  <0  =^>  kT  <  a(zz) 


Lemma  1.  [5]  a(zj )  with  odd  j  is  positive,  and  the  a(zj )  with  even  j 
is  negative. 

According  to  Lemma  1,  the  range  of  kj  becomes 

0  <  &/  <  .  min  {a(zj)}  (25) 

J=l)3,5,... 

The  above  procedure  is  summarized  in  Algorithm  1  (fully  presented 
in  the  paper  [5]). 

Applying  Algorithm  1  in  [5]  to  the  example  R  =  0.03,  N  =  80,  C  = 
3750,  we  can  get  the  stable  region  of  {kp  kp)  plane  is  showed  in  Fig.  11. 

Safe  Switching  Control  System  The  unfalsified  SSC  system,  de¬ 
scribed  in  [6],  utilizes  a  set  of  candidate  PI  controllers,  which  is  one  of 
the  most  widely  used  methods  in  the  industry.  Fundamentals  of  the  un¬ 
falsified  control  theory  and  subsequent  extension  to  the  SSC  algorithms 
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Figure  11:  Stable  region  of  (kp,kj) 


are  described  in  [6]  and  references  therein.  To  apply  this  approach,  we 
first  construct  a  bank  of  PI  controllers.  Then,  we  generate  the  ficti¬ 
tious  reference  signal  and  fictitious  error  signal  for  each  individual  PI 
controller.  Given  the  fictitious  reference  signal,  plant  input  signal,  and 
plant  output  signal  sets,  the  “best”  (optimal)  controller  is  selected  from 
the  candidate  set  using  a  properly  designed  cost  function.  To  improve 
the  overall  performance,  the  candidate  controller  parameters  can  be 
designed  off-line  using  Integral  Squared  Error  (ISE)  optimization  algo¬ 
rithms,  by  considering  the  linearized  models  in  several  general  cases. 

The  PI  controller  dynamics  can  be  expressed  as 

u  =  (kp  +  —  )e  (26) 

s 

where  e  =  y—r  is  the  queue  length  error  signal,  r  queue  length  reference, 
u  the  calculated  dropping  packet  probability,  and  y  the  actual  queue 
length. 

Cost-detectable  Cost  Function 

To  account  for  the  plant  changes  due  to  varying  N  and  R,  the  fol¬ 
lowing  ‘interval  cost  function’  is  proposed: 


Ji{t) 


max 

[tn0  it] 


e||A«||?  + 


(27) 


where  A u  is  the  deviation  of  the  dropping  probability.  tno  is  the  time 
of  the  cost  function  being  reactivated  at  the  nth  time.  The  truncated 
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time(sec) 


Figure  12:  Outputs  of  Cost-detectable  Cost  Function  (27) 

2-norm  ||x||r  =  *  ft  ( x(r))2dr . 

V  n0 

Algorithm  1  Interval  cost  function  switching  algorithm 
1:  Initialization:  define  a  set  of  candidate  controllers,  and  an  initial  controller  in  the  loop 
at  the  beginning.  Set  initial  cost  function  output  to  be  0.  Initialize  a  timer  T  =  Os. 

2:  Measure  A u  and  y.  Run  the  timer. 

3:  Calculate  r*  and  e^,  and  Ji(t). 

4:  Switch  the  controller  arg  min  JAt)  into  the  loop  if  min  JAt)+e  <  JCurrent  controller  (t). 

l<i<N  1  <i<N 

5:  Measure  y,  and  calculate  e  =  r2  —  y2. 

6:  If  |e  |  >  emax,  initialize  the  timer  to  0,  and  go  back  to  step  2. 

7:  If  |e  |  <  emax  and  T  <  tmax,  go  back  to  step  2. 

8:  If  |e  |  <  emax  and  T  >  trnax ,  stop  the  timer,  initialize  the  cost  function  output  to  be  0, 
and  go  back  to  step  5  (shut  off  the  cost  monitor). 

Remark  :  Shutting  the  cost  monitor  in  Algorithm  1  saves  the  com¬ 
putational  cost,  which  is  a  consideration  of  high  importance  in  satellite 
networks. 

Lemma  2.  [5]  Consider  the  cost  function  in  (27)  with  a  >  0.  If  the 
candidate  controllers  in  the  set  K  are  in  PI  form,  as  shown  in  (26), 
(■IK)  is  cost- detectable 

2.2.1  Simulation  Results 

The  SSC  scheme  is  simulated  using  MATLAB  and  ns-2.  The  queue 
length  reference  r  =  200  packets.  At  the  beginning,  Source  :  1  —  30 
start  sending  data  to  Destination  :  1;  at  t  =  50  s,  Source  :  1  —  60 
switch  to  send  data  to  Destination  :  1;  at  t  =  100  s,  Source  :  1  —  60 
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Figure  13:  Instant  Queue  Length  using  SSC  Scheme 


Figure  14:  Instant  Dropping  Probability  using  SSC  Scheme 


Figure  15:  Throughput  using  SSC  Scheme 
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Figure  16:  Congestion  Window  Size  using  SSC  Scheme 


switch  to  send  data  to  Destination  :  2;  at  t  =  150  s,  Source  :  1  —  120 
switch  back  to  send  data  to  Destination  :  2. 

First,  based  on  the  linearized  models  about  some  nominal  points, 
we  find  the  stable  region  using  Algorithm  1  in  [6] .  Secondly,  we  design 
the  candidate  controller  off-line,  using  ISE  optimization  algorithm,  as 
follows.  The  parameters  for  the  first  controller  are  kp  =  0.0001  and 
kj  =  0.0002  which  are  designed  around  nominal  point  N  =  60,  R  = 
0.03,  those  for  the  second  controller  are  kp  =  0.000018  and  kj  =  0.00001 
around  N  =  60,  R  =  0.26,  and  the  third  controller  is  a  fuzzy  one  as 
described  above.  In  Algorithm  1,  we  assign  emax  =  2000 packets2  and 
i'max  5  S. 

Moreover,  TCP  NewReno  scheme  is  used  in  the  detailed  simulation 
in  ns-2,  since  TCP  NewReno  is  even  closer  to  the  fluid  model  in  equa¬ 
tion  (1),  which  ignores  the  slow  start  procedure.  The  main  contribution 
of  NewReno  is  the  introduction  of  Fast-Retransmit  and  Fast-Recovery, 
which  introduces  two  indications  of  a  packet  loss:  Retransmission  Time¬ 
out  (RTO)  and  arrival  of  three  duplicate  acknowledgements.  If  a  packet 
loss  is  indicated  by  a  timeout,  TCP  believes  that  the  network  is  con¬ 
gested  and  hence  enters  into  the  slow  start  procedure  to  recover  from  it. 
Arrival  of  three  duplicate  acknowledgements  means  that  there  is  still 
data  flowing  between  the  two  ends,  and  congestion  does  not  happen. 
Hence  TCP  employs  the  Fast-Retransmit  and  Fast-Recovery  to  avoid 
reducing  the  flow  abruptly. 

We  observe  that  using  RED,  adaptive  RED,  or  single  PI  controller 
does  not  result  in  satisfactory  tracking  the  queue  length  reference.  In- 
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stead,  the  proposed  SSC  unfalsified  controller  can  track  the  reference 
signal  much  better. 

In  Fig.  12,  the  time  intervals  when  all  the  cost  function  outputs  drop 
to  0  signify  the  case  that  \e\  <  emax  and  T  >  tmax  (in  MATLAB),  so 
the  cost  monitor  is  shut  off  to  save  the  computational  cost,  and  the 
historical  data  is  discarded  to  re-initialize  the  cost  level.  The  change  of 
cost  function  outputs  from  zero  to  nonzero  signifies  the  case  that  |e  |  > 
emaX:  and  the  current  candidate  controller  might  be  falsified.  Hence,  the 
cost  monitoring  is  reactivated  to  detect  the  optimal  controller.  Fig.  14 
shows  that  the  the  dropping  probability  increases  while  the  number  of 
connections  increases  or  the  propagation  delay  decreases.  Fig.  15  shows 
that  the  throughput  in  the  whole  period  is  at  or  close  to  3750 packets/ s} 
which  means  the  bandwidth  utilization  is  pretty  high.  The  congestion 
window  size  of  one  TCP  flow  is  zoomed  in  and  showed  in  Fig.  16,  which 
is  consistent  with  the  TCP  NewReno  scheme. 
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2.3  Cross  Entropy  Accelerated  Ant  Routing  System  for  Nano- 
Satellite  Networks 

The  previous  section  derived  some  new  congestion  control  algorithms; 
this  section  develops  SSA  data  routing  techniques.  An  efficient  routing 
algorithm  plays  a  key  role  in  optimizing  network  resources  in  a  LEO 
satellite  network.  We  have  developed  a  new  routing  system  based  on 
the  Cross  Entropy  optimization  method  (a  heuristic  method  for  solv¬ 
ing  combinatorial  optimization  problems)  [3] ,  [4] .  The  novel  on-demand 
routing  method,  named  Cross  Entropy  Accelerated  Ant  Routing  Sys¬ 
tem  (CEAARS)  is  applied  in  simulations  to  the  regular  constellation 
LEO  satellite  networks.  By  implementing  simulations  on  an  Iridium- 
like  satellite  network,  we  compared  the  proposed  CEAARS  algorithm 
with  two  standard  approaches  to  adaptive  routing  protocols  on  the  In¬ 
ternet:  Distance- Vector  (DV)  and  Link-State  (LS),  as  well  as  with  the 
original  Cross  Entropy  Ant  Routing  System  (CEARS).  DV  algorithms 
are  based  on  the  distributed  Bellman  Ford  algorithm,  and  LS  algo¬ 
rithms  are  an  implementation  of  Dijkstra’s  single  source  shortest  path 
algorithm.  The  results  show  that  CEAARS  not  only  remarkably  im¬ 
proves  the  convergence  speed  of  achieving  optimal  or  suboptimal  paths, 
but  also  reduces  the  number  of  overhead  ants  (management  packets). 

Since  the  distance  between  satellite  planes  changes  with  the  move¬ 
ment  of  the  satellites,  the  interplane  intersatellite  links  (ISLs)  are  longest 
when  satellites  are  over  the  equator  and  shortest  when  they  are  over 
the  polar  region  boundaries.  In  the  near  polar  (or  Walker  star)  constel¬ 
lations,  such  as  Iridium  and  Teledesic,  interplane  ISLs  are  not  main¬ 
tained  as  satellites  come  close  to  the  poles,  due  to  adverse  pointing 
and  tracking  conditions,  and  are  reestablished  when  satellites  move 
to  lower  latitudes.  Since  the  satellite  movements  cause  changes  in  the 
network  topology,  more  robust  routing  algorithms  for  satellite  networks 
are  needed,  which  can  establish  and  maintain  an  optimal  path  between 
source  and  destination  in  a  dynamically  varying  topological  setting. 

Prior  results  in  the  routing  algorithms  for  satellite  networks  mainly 
focus  on  the  static  snapshot  path  setup  and  snapshot  transition  phases. 
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Satellites  forward  the  packets  according  to  their  routing  tables,  which 
are  generated  in  a  ground  switch  centrally,  and  the  routing  tables  are 
updated  according  to  a  switching  table,  as  the  satellite  network  topol¬ 
ogy  changes.  In  the  above  schemes,  adaptivity  and  robustness  of  a 
satellite  networks  are  not  optimal.  Moreover,  high  energy  consump¬ 
tion  is  required  to  maintain  all  snapshot  routes  subject  to  the  variable 
topology  requirements,  which  is  inappropriate  for  satellite  networks  (es¬ 
pecially  nano-satellite  networks  with  strict  energy  constraint). 

In  the  present  work,  a  robust  Cross  Entropy  Ant  Routing  Sys- 
tem(CEARS)  and  its  improved  version  Cross  Entropy  Accelerated  Ant 
Routing  System(CEAARS)  are  introduced  for  regular  constellations 
LEO  satellite  networks  (such  as  near  polar  -or  Walker  star-  constel¬ 
lations),  and  compared  with  the  distance- vector  (DV)  and  link-state 
(LS),  which  are  two  known  approaches  to  adaptive  routing  protocols 
on  the  Internet  Protocol  (IP)  networks.  Both  CEARS  and  CEAARS 
are  adaptive  to  the  change  of  network  topology,  and  also  robust  to  the 
path  recovery  due  to  the  link  and  node  failure.  CEARS  is  based  on  the 
global  random  search,  so  it  takes  a  much  longer  time  to  find  an  opti¬ 
mal  or  suboptimal  path  at  the  beginning,  and  consumes  more  energy. 
Due  to  the  capability  of  sensing  the  direction  of  destination,  CEAARS 
is  characterized  by  a  faster  convergence  speed  and  better  energy  con¬ 
sumption. 

An  Iridium- like  constellation,  shown  in  Fig.  17,  which  consists  of  66 
satellites  with  6  planes  of  11  satellites  per  plane  in  a  near  polar  LEO,  is 
considered  as  an  example.  Assume  that  each  satellite  or  node  is  capable 
of  having  four  links,  comprising  two  intraplane  ISLs  forward  and  back 
and  two  interplane  ISLs  to  the  satellites  in  the  left  and  right  neighboring 
planes.  However,  cross  seam  ISLs  are  inactive,  and  interplane  ISLs  are 
inactive  in  the  near-polar  region.  This  topology  can  be  considered  as 
a  variation  of  the  bidirectional  seamless  Manhattan  network,  as  shown 
in  Fig.  18. 

By  examining  the  number  of  paths,  we  show  in  [3]  that  it  is  a  rare 
event  to  find  an  optimal  path  by  doing  a  random  walk. 

In  the  first  phase  of  the  Cross  Entropy  Ant  Routing  System  (CEARS) 
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Figure  17:  Layout  of  Iridium  constellation 


method,  ants  (or  agents)  walk  randomly  through  the  network  from  the 
source  node  to  the  destination  node,  measure  the  cost  of  the  paths  and 
deposit  pheromone  to  represent  the  quality  of  the  traversed  paths.  In 
the  second  phase,  according  to  the  cost  of  the  discovered  paths,  the 
cross  entropy  between  the  density  of  the  cost  of  paths  and  the  density 
of  the  routing  probabilities  is  minimized  to  generate  the  optimal  (lowest 
cost)  path,  and  the  routing  probability  matrix  is  updated  accordingly. 


2.3.1  Theory 

‘Close  to  Destination’  Neighbor  Nodes  Set. 

The  following  simulation  results  show  that  CEARS  (which  does  not 
take  into  account  the  relationship  among  the  locations  of  the  satellite 
nodes)  takes  a  relatively  long  time  to  converge  to  the  optimal  paths. 
This  makes  the  CEARS  not  suitable  as  the  routing  algorithm  for  short- 
time  connections.  Considering  the  fact  that  the  satellite  nodes  have 
certain  relation  with  each  other  in  regular  constellations,  we  want  to 
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Figure  18:  Manhattan  Network  Topology 
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improve  the  convergence  speed  by  adding  the  so-called  destination- 
oriented  ants,  whose  neighbor  set  includes  only  close-to-destination 
neighbor  nodes. 

The  close-to-destination  neighbor  nodes  set  is  a  set  in  which  the 
nodes  are  contained  in  the  minimum  hop  path.  To  determine  the  close- 
to-destination  neighbor  nodes  set  with  minimum  hop  metrics,  we  use 
the  location  relations  between  the  current  satellite  nc  and  the  des¬ 
tination  satellite  nr,  as  shown  in  the  pseudo  code  in  [4],  in  which  the 
neighbor  nodes  set  is  narrowed  down  by  comparing  the  row  and  column 

indices  of  the  current  satellite  with  those  of  the  destination  satellite. 

The  CEAARS  algorithm  is  described  as  follows. 

1.  A  uniformly  distributed  probability  matrix  is  initialized.  Three  types  of  ants:  destination- 
oriented  ants,  normal  ants  and  explorer  ants,  are  generated  at  the  source  node  and 
start  to  search  forward  paths  to  the  destination  node. 

2.  For  the  destination-oriented  ants,  according  to  the  relative  position  between  the  cur¬ 
rent  node  and  destination  node  by  comparing  their  row  and  column  indexes,  a  set 
of  “close-to-destination”  neighbor  nodes  are  generated.  For  other  types  of  ants,  the 
neighbor  set  includes  all  neighbor  nodes. 

3.  For  normal  ants  and  destination-oriented  ants,  the  next  hop  node  to  be  visited  is 
selected  according  to  a  probability  distribution  in  their  own  neighbor  sets.  If  no  nodes 
were  available  in  the  set  of  “close-to-destination”  neighbor  nodes  for  destination- 
oriented  ants,  the  set  with  all  neighbor  nodes  will  be  considered.  For  explorer  ants, 
the  next  hop  node  to  be  visited  is  selected  randomly. 

4.  If  the  selected  node  has  been  visited,  it  will  be  avoided. 

5.  If  the  forward  ants  does  not  arrive  at  the  destination,  go  to  (2). 

6.  If  the  forward  ants  arrived  at  the  destination,  the  cost  of  the  path  is  calculated,  and 
the  temperature  is  updated. 

7.  Backward  ants  are  generated  and  return  on  the  reserved  path,  updating  the  pheromone 
values  and  probability  distribution  values. 

Fuzzy  Adaptive  Ant  Generation  Rates  In  CEARS,  a  relatively 
high  ant  generation  rate  is  necessary  to  maintain  up-to-date  statistics  of 
the  network  status.  This  might  create  significant  overhead  in  terms  of 
the  routing  traffic,  increasing  the  cost  and  having  a  negative  impact  on 
the  overall  network  performance.  To  save  computation  and  bandwidth 
cost  in  a  satellite  network,  ant  generation  rates  are  adaptively  changed 
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according  to  the  routing  convergence  condition.  For  example,  at  the 
beginning,  all  forward  ants  visit  the  next  hop  node  randomly,  and  the 
backward  ants  follow  their  own  forward  paths  reversely  and  update 
the  pheromone  values.  During  this  period,  the  backward  rate  is  much 
lower  than  the  forward  rate.  In  this  case,  we  want  to  assign  a  high  ant 
generation  rate  to  find  the  optimal  paths  as  soon  as  possible.  Then, 
more  ants  follow  the  lowest  cost  paths  forward  and  backward,  and  the 
backward  rate  becomes  close  to  the  forward  rate.  In  this  case,  we  want 
to  assign  a  low  ant  generation  rate  to  save  computation  and  bandwidth 
cost.  Even  though  we  know  how  to  control  the  ant  generation  rate  by 
using  human  heuristic  knowledge,  it  is  difficult  to  use  formal  analysis  to 
calculate  the  explicit  ant  generation  rate.  Hence  we  adopt  fuzzy  logic, 
which  is  easy  to  use  and  very  flexible. 

For  the  convenience  of  calculation  and  satisfactory  accuracy,  we  di¬ 
vide  A C,  rng ,  rS5,  regi  and  (3  into  2  fuzzy  sets:  Big  (B)  and  Small 
(S),  where  rngjrsg  and  reg  represent  normal  ants,  destination-oriented 
ants  and  explorer  ants  generation  rates  respectively.  AC  is  the  dif¬ 
ference  between  the  average  cost  and  a  pre-set  reference  cost.  In  this 
work,  the  non-optimal  actual  cost  is  usually  much  greater  than  the 
theoretical  minimum  cost  (more  than  ten  times  greater).  Using  the 
fuzzy  logic  rules,  we  can  define  the  reference  cost  to  be  three  times 
the  minimum  cost,  and  consider  all  lower  cost  as  ’SMALL’,  otherwise 
’BIG’.  By  rf{k)  and  rb(k)  we  denote  forward  rate  and  backward  rate 
at  visit  k,  respectively,  and  rbf  =  Then  we  can  use  the  heuristic 
knowledge  about  the  relationship  between  forward  rate  and  backward 
rate  to  design  some  simple  fuzzy  rules  for  the  ant  generation  rates  and 
evaporation  rate  shown  in  Tables  1  and  2. 

The  control  rules  are  developed  with  ri,f  and  AC  as  premise  vari¬ 
ables,  and  rngi  rsg ,  reg ,  and  (3  as  consequent  variables  of  the  each  rule. 
The  structure  of  the  fuzzy  rule  is  described  as: 


IF  rbf  is  SMALL  and  AC  is  BIG, 

THEN  rng  and  rsg  are  BIG,  reg  is  SMALL,  (3  is  BIG. 
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In  this  work,  the  corresponding  fuzzy  sets  for  r&f,  A C,  rng,  rsg ,  reg  and 
/ 3  are  specihcally  {0.3,  0.6},  {0.05,0.1},  {40,100},  {40,100},  {5,20}, 
and  {0.6,0.95}.  The  normal  triangular  membership  function  is  used. 
For  the  defuzzification  procedure  that  produces  a  crisp  control  output 
from  the  result  of  the  inference,  the  Center  of  Gravity  method  is  used 
here  to  get  the  final  output  value. 


Table  1:  Rule  base  for  ant  generation  rates 


reg 

rbf 

S 

B 

AC 

s 

s 

s 

B 

s 

B 

rng  or  rsg 

rbf 

S 

B 

AC 

S 

S 

S 

B 

B 

S 

In  order  to  play  down  the  impact  of  the  bursty  traffic  on  the  forward 
and  backward  ants  rates,  an  exponentially  weighted  moving  average 
(EWMA)  low  pass  filter  is  designed  to  estimate  the  rates  as: 

rf{k)  <-  Wrrf(k  -  1)  +  (1  -  Wr)rg{k  -  1) 
n(k) «-  Wrrb(k  - 1)  +  (l  -  wr)— 

T 

where  Wr  is  the  time  constant  for  the  low  pass  filter,  r&( 0)  =  0,  and 
is  the  number  of  backward  ants  in  time  [(k  —  1  )r,  kr\. 

2.3.2  Simulation  Results 

For  performance  evaluation  of  the  CEAARS,  simulation  experiments 
on  an  Iridium-like  topology  network  are  studied  on  the  Network  Sim¬ 
ulator  2  (ns-2),  a  discrete  event  simulator.  The  following  assumptions 
are  considered:  the  propagation  delay  of  each  ISL  is  15ms,  and  all  in¬ 
terplane  ISLs  in  one  interplane  become  inactive  simultaneously  in  the 
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Table  3:  Interplane  ISLs  status  on  different  Snapshots 


Si 

*5*2 

S3 

s4 

S5 

So 

Ay 

Ag 

A9 

A,0 

An 

h 

0 
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End  to  End  Delay 


Figure  19:  Average  end  to  end  delay  using  DV  (Scenario  1) 
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Average  end  to  end  delay  (s) 


End  to  End  Delay 


Figure  20:  Average  end  to  end  delay  using  LS  (Scenario  1) 
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End  to  End  Delay 


Figure  21:  Average  end  to  end  delay  using  GEARS  (Scenario  1) 
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End  to  End  Delay 


Figure  22:  Average  end  to  end  delay  using  CEAARS  (Scenario  1) 
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Table  4:  Quality  of  paths  found  in  simulation  Scenario  1  Pattern  1  at  different 
snapshots  using  the  three  different  schemes.  Values  are  obtained  from  the  simulation 
after  the  first  11  snapshots  (550  seconds)  in  the  form  of  Mean(Standard  Deviation) 
seconds 


Snap 

Short  Delay 

CEARS  Delay 

CEAARS  Delay 

DV  Delay 

LS  Delay 

Si 

0.0750 

0.1208(0.0128) 

0.0905(0.0000) 

0.3541(0.0404) 

0.4647(0.0754) 

s2 

0.0750 

0.1430(0.0561) 

0.0899(0.0009) 

0.3692(0.0023) 

0.5010(0.0037) 

S3 

0.0750 

0.1292(0.0199) 

0.0904(0.0000) 

0.3743(0.0024) 

0.5028(0.0029) 

S4 

0.1050 

0.1625(0.0197) 

0.1302(0.0059) 

0.3952(0.0414) 

0.5296(0.0410) 

S5 

0.1050 

0.1557(0.0090) 

0.1308(0.0020) 

0.4108(0.0027) 

0.5314(0.0070) 

S(j 

0.1350 

0.2296(0.0611) 

0.1649(0.0054) 

0.4351(0.0419) 

0.4336(0.0407) 

S7 

0.1050 

0.2057(0.0204) 

0.1295(0.0071) 

0.4503(0.0028) 

0.4524(0.0040) 

S8 

0.1050 

0.1728(0.0379) 

0.1263(0.0008) 

0.4577(0.0025) 

0.4553(0.0041) 

S9 

0.0750 

0.1314(0.0274) 

0.0928(0.0065) 

0.3542(0.0447) 

0.4592(0.0172) 

Sio 

0.0750 

0.1206(0.0092) 

0.0905(0.0000) 

0.3702(0.0023) 

0.4654(0.0040) 

Sn 

0.0750 

0.1246(0.0179) 

0.0905(0.0000) 

0.3750(0.0023) 

0.4668(0.0028) 

Table  5:  Quality  of  paths  found  in  simulation  Scenario  1  Pattern  4  at  different 
snapshots  using  the  three  different  schemes.  Values  are  obtained  from  the  simulation 
after  the  first  11  snapshots  (550  seconds)  in  the  form  of  Mean(Standard  Deviation) 
seconds 


Snap 

Short  Delay 

CEARS  Delay 

CEAARS  Delay 

DV  Delay 

LS  Delay 

Si 

0.0750 

0.0909(0.0056) 

0.0764(0.0000) 

0.0777(0.0000) 

0.0777(0.0000) 

s2 

0.0750 

0.1187(0.0607) 

0.0763(0.0000) 

0.0777(0.0000) 

0.0777(0.0000) 

S3 

0.0750 

0.0932(0.0062) 

0.0764(0.0000) 

0.0777(0.0000) 

0.0777(0.0000) 

S4 

0.1050 

0.1259(0.0233) 

0.1108(0.0051) 

0.1087(0.0003) 

0.1088(0.0016) 

S5 

0.1050 

0.1219(0.0068) 

0.1121(0.0020) 

0.1087(0.0000) 

0.1087(0.0006) 

S6 

0.1350 

0.1535(0.0134) 

0.1391(0.0045) 

0.1398(0.0003) 

0.1398(0.0000) 

Sr 

0.1050 

0.1295(0.0141) 

0.1105(0.0059) 

0.1093(0.0043) 

0.1093(0.0042) 

s8 

0.1050 

0.1183(0.0054) 

0.1082(0.0007) 

0.1087(0.0000) 

0.1087(0.0000) 

s9 

0.0750 

0.1034(0.0140) 

0.0781(0.0049) 

0.0783(0.0043) 

0.0783(0.0043) 

Sio 

0.0750 

0.0926(0.0065) 

0.0764(0.0000) 

0.0777(0.0000) 

0.0777(0.0000) 

Sn 

0.0750 

0.0907(0.0057) 

0.0764(0.0000) 

0.0777(0.0000) 

0.0777(0.0000) 
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Table  6:  Total  routing  overhead  packets  (packets)  under  Scenario  1 


Pattern  1 

Pattern  2 

Pattern  3 

Pattern  4 

DV 

73874 

74036 

74034 

74046 

LS 

343120 

277266 

342416 

466161 

CEARS 

66085 

66101 

66194 

66105 

CEAARS 

36649 

37708 

36758 

36232 

polar  region.  To  clearly  show  the  simulation  results,  we  divide  the  net¬ 
work  into  11  snapshots  based  on  the  status  of  interplane  ISLs,  and  one 
snapshot  lasts  50  seconds.  Table  3  shows  the  interplane  ISL’s  status, 
where  Si  denotes  the  ith  snapshot;  Ij  denotes  the  jth  interplane;  0  at 
the  i  —  th  column  and  jth  row  means  that  the  interplane  ISLs  of  the  jth 
interplane  at  the  ith  snapshot  are  inactive,  and  1  means  active.  In  the 
simulation  the  source  node  is  3,  and  the  destination  node  is  58.  Maxi¬ 
mum  two  hundred  data  packets  with  the  size  of  1000  bytes  per  packet 
are  generated  per  second  at  the  source  node.  For  CEARS,  (5  =  0.95, 
and  100  normal  ants  and  10  explorer  ants  are  generated  per  second  at 
the  source  node. 

With  the  configuration  explained  above,  CEAARS  has  been  imple¬ 
mented  and  compared  with  the  Distance  Vector(DV)  and  Link  State 
(LS)  approaches,  under  both  light  and  heavy  load  conditions  as  follows. 

Scenario  1  ( Relation  between  load  and  performance ): 

—  Pattern  1  ( Heavy  load):  the  ISL  bandwidth  is  set  as  1.5 Mbit/s. 

—  Pattern  2:  the  ISL  bandwidth  is  set  as  5Mbit/ s. 

—  Pattern  3:  the  ISL  bandwidth  is  set  as  10 Mbit/s. 

—  Pattern  4  ( Light  load):  the  ISL  bandwidth  is  set  as  15Mbit/ s. 

By  inspection  of  the  topology,  we  can  calculate  the  minimum  prop¬ 
agation  delay  from  the  source  node  to  the  destination  node  at  the  first 
and  last  three  snapshots  as  75ms,  at  the  4 th,  5th,  7th  and  8th  snap¬ 
shots  as  105ms,  and  at  the  6 th  snapshot  as  135ms.  Figures  19,  20,  21 
and  22  show  the  average  end  to  end  delay  results  in  Scenario  1.  We 
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can  find  that,  at  the  first  50  seconds,  the  resulting  delay  using  CEARS 
is  around  400ms  which  is  quite  high  corresponding  to  the  theoretical 
minimum  delay  at  the  two  snapshots,  while  DV,  LS,  and  CEAARS 
converge  to  the  optimal  paths  very  fast  when  the  traffic  load  is  light. 
DV  and  LS  always  select  one  of  the  shortest  paths  to  forward  packets. 
In  the  light  load  traffic  situation,  selecting  one  of  the  shortest  paths  to 
forward  packets  is  the  best  selection.  However,  in  the  heavy  load  traffic 
situation,  DV  and  LS  cannot  avoid  the  congested  path.  Table  4  and 
Table  5  show  the  statistical  measures  of  the  mean  and  standard  devi¬ 
ation  values  in  the  extended  simulation.  By  comparing  Figs.  19,  20, 
and  22,  also  from  Table  4,  we  can  see  that  CEAARS  outperforms  DV 
and  LS  in  heavy  load  situation,  since  CEAARS  can  adaptively  select 
the  minimum  cost  path.  Moreover,  from  Table  6,  we  can  find  that 
the  number  of  routing  overhead  packets  in  LS  is  much  higher  than 
in  DV,  CEARS  and  CEAARS,  because  the  new  state  information  in 
LS  needs  to  be  disseminated  throughout  the  network  by  the  LS  flood¬ 
ing  procedure.  In  DV  algorithms,  each  node  uses  neighbors’  routing 
tables  to  get  the  information  of  link  costs  to  all  its  neighbors,  and 
trigger  update  starts  only  when  link  changes  are  affecting  the  shortest 
path.  This  saves  a  lot  of  overhead  packets  especially  in  LEO  satellite 
networks,  where  the  ISLs  become  alternately  inactive  and  active.  In 
CEARS  and  CEAARS,  backward  ants  update  the  routing  table  at  in¬ 
tervals  and  there  is  no  obligation  to  have  global  updating  mechanisms. 
Overall,  CEAARS  demonstrates  good  adaptiveness  to  topology  vari¬ 
ations  in  simulations.  Unlike  DV  and  LS,  CEAARS  can  adaptively 
choose  the  optimal  path  based  on  a  cost  function  with  no  obligation  to 
global  updating  mechanisms.  Since  CEAARS  is  capable  of  sensing  the 
direction  of  destination,  it  remarkably  improves  the  convergence  speed. 
Since  CEAARS  can  find  the  optimal  or  suboptimal  paths  quite  soon 
at  the  beginning  of  the  simulation,  it  is  suitable  for  both  short-time 
and  long-time  connections.  The  simulation  on  CEAARS  also  shows  a 
significant  saving  in  the  number  of  ants  with  similar  data  packet  delay 
and  service  availability.  Consequently,  CEAARS  can  be  implemented 
as  an  on-demand  routing  system  with  low  memory  and  computations. 
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2.4  Satellite  Network  TCP  based  Congestion  Control  using 
Convex  Optimization 


In  addition  to  transmitting  the  SSA  data  across  a  single  layer  (such  as 
LEO),  the  data  also  must  often  be  transmitted  between  layers  (such  as 
to  GEO).  A  multi-layered  satellite  network  consisting  of  a  large  number 
of  nano-satellites  is  ideally  suited  to  perform  space  situational  aware¬ 
ness.  The  large  number  of  nano  satellites  collects  visual  information  of 
space  objects  and  transfer  data  to  ground  stations  through  a  layered 
network  in  real  time.  Due  to  the  dynamic  topology  of  the  network, 
large  propagation  delays  and  bulk  data  transfers,  significant  delays  and 
loss  of  data  can  be  encountered  with  traditional  TCP  congestion  control 
schemes.  To  overcome  these  drawbacks,  network  snapshots  consisting 
of  layered  dynamic  clusters  of  satellites  are  formed  with  modified  TCP 
congestion  control  schemes  such  that  the  delay  and  packet  losses  are  at 
a  minimum.  The  layered  dynamic  clusters  are  allowed  to  transmit  data 
at  maximum  possible  rates  till  congestion  is  detected.  On  congestion, 
an  optimized  data  rate  is  determined  using  Convex  Optimization  tech¬ 
niques  to  support  high  data  rates  and  reduce  congestion.  The  multi¬ 
layered  satellite  architecture  comprises  of  three  layers.  The  lowest  layer 
comprises  of  nano-satellites  (N)  with  limited  communication  capabili¬ 
ties.  Layer  2  comprises  of  MEO  satellites  (M)  and  layer  3  comprises 
of  GEO-stationary  satellites  (G).  Satellites  communicate  within  layers 
via  links  known  as  the  intra-orbital  links  whereas  satellite  communica¬ 
tions  between  layers  is  done  using  the  links  known  as  inter-orbital  links. 
Inter-orbital  data  flow  is  assumed  to  exist  only  in  the  lower  and  GEO 
layer.  Only  the  GEO  stationary  satellites  are  assumed  have  the  capa¬ 
bilities  to  communicate  with  the  ground  stations.  The  multi-layered 
satellite  network  is  assumed  so  that  the  top  layer  will  always  represent 
the  GEO  layer  and  the  other  two  layers  can  be  located  anywhere  from 
the  low  earth  orbits  to  right  below  GEO  satellites.  Next,  the  multi¬ 
layered  architecture  is  further  modified  into  a  layered  dynamic  cluster 
of  satellites  by  selecting  network  snapshots  manager  or  cluster  head 
satellites  in  each  layer.  First,  several  nano-satellites  in  the  lowest  layer 
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are  chosen  as  the  dynamic  cluster  heads  based  on  certain  communi¬ 
cation  capabilities  and  locations  of  neighboring  nano-satellites.  These 
cluster  heads  will  collect  data  from  its  members  in  the  lower  layer  and 
forward  the  data  to  a  single  satellite  in  the  layer-2.  Each  satellite  in  the 
layer- 2  will  assume  a  role  of  dynamic  cluster  head  however  will  commu¬ 
nicate  only  with  dynamic  cluster  heads  in  the  lower  layer  and  forward 
the  received  data  to  an  appropriate  GEO  layer  satellite.  At  the  dy¬ 
namic  cluster  heads  in  each  layer  the  data  is  accumulated.  Depending 
on  the  incoming  data  rate,  size  of  the  data  buffers  and  the  outgoing 
data  rate,  significant  queuing  delay  would  be  introduced  leading  to  drop 
in  packets.  A  simple  solution  to  avoid  packet  drops  would  be  to  have 
a  large  data  buffer  which  is  not  practically  feasible.  Even  if  sufficient 
data  buffers  are  available  the  queuing  delay  would  still  be  significant 
due  to  low  outgoing  data  rates.  The  other  simple  solution  is  to  re¬ 
duce  significantly  the  incoming  data  rates  which  in  turn  would  increase 
the  end  to  end  delay  in  transmitting  the  data  to  the  ground  station. 
Hence,  there  is  a  need  for  determining  an  optimal  layered  cluster  and 
optimum  incoming  data  rates.  In  this  work,  using  convex  optimization 
techniques,  optimal  incoming  data  rates  are  computed  at  each  satellite 
irrespective  of  their  layers.  They  are  computed  in  each  snapshot  to 
keep  the  delay  and  packet  loss  at  a  minimum.  During  the  process  of 
determining  optimal  incoming  data  rates  a  collection  of  satellites  are 
also  identified  to  perform  a  particular  space  situational  awareness  task. 
In  the  next  section,  some  of  the  results  are  presented  and  detailed  dis¬ 
cussions  on  the  TCP  modeling  with  convex  optimization  are  provided 
in  the  appendix. 

2.4.1  Simulation  Results 

A  network  comprising  of  120  Nano-satellites  in  the  lowest  layer,  12 
satellites  in  the  middle  layer  and  3  GEO  satellites  were  considered.  It 
was  assumed  that  only  one  of  these  GEO  satellites  has  connectivity  to 
the  ground  station  located  at  Florida.  The  other  two  GEO  satellites 
forward  their  data  to  the  GEO  satellite  connected  to  the  ground.  This 
scenario  was  selected  to  stress  the  system  and  thereby  demonstrate  the 
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effect  on  buffer  usage,  bandwidth  utilization,  etc.  The  nano-satellites 
in  the  lower  layer  are  assumed  to  have  a  transmission  capability  of  one 
tenth  that  of  the  satellites  in  the  middle  layer.  The  GEO  satellites 
are  assumed  to  have  the  highest  transmission  capability.  However,  the 
physical  distance  (and  in  turn  the  propagation  delay)  between  lower  to 
middle  layers  and  middle  to  GEO  layers  are  significantly  smaller  than 
the  propagation  delays  between  GEO  satellites  and  GEO  satellite  to 
the  ground  station.  In  a  snapshot  a  GEO  satellite  could  communicate 
with  a  maximum  of  four  satellites  in  the  middle  layer  and  a  satellite  in 
the  middle  layer  could  communicate  with  a  maximum  of  three  cluster 
head  satellites  in  the  lower  layer.  The  cluster  size  in  the  lower  layer  was 
restricted  to  five  including  the  cluster  head.  Simulations  were  run  for  a 
fixed  amount  of  time  and  performance  parameters  at  each  transmission 
were  collected.  The  simulation  run  included  a  number  of  snapshots 
with  fixed  snapshot  duration.  It  is  assumed  that  there  is  no  change  in 
positions  of  the  satellite  during  a  snapshot.  After  a  snapshot,  satellites 
recalculate  their  positions,  memberships/leadership  and  data  rates.  In 
Figures  23,  24,  and  25,  the  average  of  the  sum  of  buffer  utilization  of 
all  cluster  head  satellites  for  each  of  the  layers  with  snapshot  duration 
of  120  secs  and  with  traditional  TCP  congestion  control  are  shown. 

From  Figure  23,  it  can  be  seen  that  the  buffer  utilization  of  the 
cluster  head  nano-satellites  is  around  20%.  This  is  due  to  the  small 
propagation  delay  between  the  lower  and  middle  layer.  The  buffer  uti¬ 
lization  of  the  satellites  in  the  middle  layer  is  negligible  due  to  the 
higher  transmission  rate  capability  between  the  middle  layer  and  the 
GEO  layer.  However,  in  case  of  the  GEO  satellites  as  shown  in  Figure 
25,  the  buffer  utilization  is  greater  than  100%.  This  is  due  to  the  large 
propagation  delay  between  the  GEO  satellites  and  the  ground  station. 
This  will  introduce  a  significant  delay  and  cause  loss  of  packets.  In 
Figure  26,  the  average  of  the  sum  of  bandwidth  utilization  at  GEO 
satellites  is  shown.  In  Figure  26,  it  can  be  seen  that  the  bandwidth 
utilization  consistently  falls  below  100%  due  to  the  TCP  congestion 
control  trying  to  adjust  the  window  size  through  its  built  in  slow  start 
and  congestion  control  mechanisms.  In  Figure  27,  the  average  of  the 
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25 


Figure  23:  Average  of  Sum  of  Buffer  Utilization  of  all  cluster  head  nano-satellites 
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Figure  24:  Average  of  Sum  of  Buffer  Utilization  of  all  satellites  in  the  middle  layer 
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Average  of  Sum  of  Buffer  Utilization  of  all  Geo  Satellites  In  % 


Figure  25:  Average  of  Sum  of  Buffer  Utilization  of  all  satellites  in  the  GEO  layer 
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Average  of  Sum  of  Bandwidth  Utilization  of  all  Geo  Satellites  In 


Number  cf  Transmissions 


Figure  26:  Average  of  Sum  of  Bandwidth  Utilization  of  GEO  satellites 
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Figure  27:  Average  of  Sum  of  Buffer  Utilization  of  all  satellites  in  the  GEO  layer 


sum  of  buffer  utilization  of  all  GEO  satellites  with  snapshot  duration 
of  120  secs  and  with  convex  optimization  based  TCP  congestion 
control  is  shown.  Comparing  Figures  25  and  27,  it  can  be  seen  that  the 
buffer  utilization  is  always  around  50%  with  the  convex  optimization 
based  TCP  congestion  control.  The  spikes  in  buffer  utilization  are  due 
to  the  satellites  recalculating  their  positions,  memberships/leadership 
and  data  rates  at  the  beginning  of  a  new  snapshot.  In  a  new  snapshot 
at  the  beginning  the  traditional  TCP  congestion  control  is  used.  Once 
the  congestion  reaches  a  critical  level,  the  convex  optimization  based 
control  reduces  the  congestion.  This  can  be  seen  in  Figure  28,  which 
shows  the  convex  optimization  mechanism  operating  in  a  single  snap¬ 
shot.  In  Figure  28,  it  can  be  seen  that  the  buffer  utilization  is  high 
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Average  of  Sum  of  Buffer  Utilization  of  all  Nano-eatellltee  In  % 


Figure  28:  Buffer  Utilization  of  GEO  Satellites  in  a  120  second  Snapshot 
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Figure  29:  Average  of  Sum  of  Bandwidth  Utilization  of  GEO  satellites  with  Convex 
Optimization  based  TCP  Congestion  Control 

at  the  beginning  of  the  snapshot  due  to  the  traditional  TCP  conges¬ 
tion  control.  After  the  congestion  reaches  the  critical  level,  the  convex 
optimization  based  TCP  congestion  control  kicks  in  and  maintains  the 
buffer  utilization  level  around  50%  for  the  entire  duration  of  the  snap¬ 
shot.  In  Figure  29,  the  average  of  the  sum  of  bandwidth  utilization  at 
GEO  satellites  with  convex  optimization  based  TCP  congestion  con¬ 
trol  is  shown.  Comparing  Figures  27  and  29,  it  can  be  noticed  that  the 
GEO  satellites  with  convex  optimization  based  TCP  congestion  control 
ensures  a  constant  100  %  bandwidth  usage  during  a  snapshot. 
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2.5  Models  for  Reconstruction  of  Images  from  Multiple  Views 


While  most  of  this  study  concentrated  on  how  SSA  data  can  be  col¬ 
lected  and  transmitted,  we  also  studied  basic  signal  processing  methods 
suited  for  this  SSA  task.  Reconstruction  of  detailed  satellite  images 
from  multiple  views  and  assembly  of  these  views  from  orbital  observa¬ 
tions  require  one  to  ask:  Under  which  transformations  do  band-limited 
signals/images  remain  invariant?  This  question  was  fully  answered  in 
our  study  sponsored  by  the  AFOSR  in  [7].  In  this  paper,  it  is  shown 
that  these  spaces  remain  invariant  if  and  only  if  the  transformations  are 
affine.  Equivalently  stated,  the  class  of  warps  preserving  Paley-Wiener 
spaces  are  limited  to  translations  and  dilations.  To  account  for  aberra¬ 
tions  and  nonlinear  warps,  it  was  subsequently  asked:  What  happens  if 
the  signal  deformations  are  not  affine?  What  are  the  function  spaces  to 
which  the  range  signals/images  belong  and  can  one  provide  an  inversion 
formula  to  recover  the  pre-images  of  these  maps?  To  address  this  prob¬ 
lem,  structure  theorems  were  constructed  to  realize  images  as  members 
of  certain  reproducing  kernel  Hilbert  spaces,  where  the  reproducing 
kernels  account  for  the  aberration.  In  this  orthogonal  polynomial  set¬ 
ting,  best  approximation  to  the  signal/image  is  possible  and  an  explicit 
inversion  formula  is  derived.  These  results  were  also  generalized  from 
the  Paley-Wiener  spaces  to  the  deBranges-Rovnyak  spaces,  which  are 
the  weighted  analogues  of  the  former  spaces.  The  weights  arise  where 
illumination  is  not  ideal  and  various  obstructions  are  present.  State¬ 
ment  of  the  exact  theorems  in  [7]  would  introduce  significant  notational 
complexity.  The  referenced  paper  is  intended  to  serve  as  a  susbsitute. 

Having  determined  the  invariance  properties  of  the  image  spaces  and 
investigated  best  approximation  approaches  to  dealing  with  the  defor¬ 
mation  and  warping  of  these  images,  the  question  of  feature  tracking 
and  realization  of  the  warp(s)  became  essential.  To  describe  the  mod¬ 
eling,  we  shall  use  /  :  M3  — >•  M2x2  to  denote  the  view  of  a  3-D  RSO, 
Q  C  1J,  by  a  camera.  Depending  on  the  model  used  for  the  camera, 
and  whether  optical,  laser,  RF  or  other  methods  are  used,  /  can  be 
appropriately  adapted.  If  h  :  D  x  [0,  oc)  — >  M3  is  a  time-dependent 
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diffeomorphism  of  D  into  M3  then  a  model  for  the  deformation  between 
two  images  of  the  same  scene  is  given  by 

Ih(x,t)  =  I{h{x,t ))  Vxeflx  n  h~\n)  c  M2x2.  (28) 

The  image  Ih{x,t )  is  the  evolving  image  of  the  3-D  RSO,  moving  at 
a  relatively  high  speed  relative  to  a  camera,  as  the  RSO  rotates  and 
translates  through  space.  If  R(t)  :  [0  G  oo)  G  50(3),  and  T(t)  : 
[0,  oo)  -G  M3  relative  to  the  Hill’s  frame,  then  h(x,t )  =  R(t)x  +  T(t) 
describes  an  affine  deformation  of  the  3-D  object.  In  this  affine  model, 
points  in  each  image  do  not  undergo  the  same  motion  but  the  motion 
of  each  point  is  related  to  its  original  location  via  an  affine  map.  This 
model  is  a  good  approximation  for  small  planar  sections  parallel  to  the 
image  plane  undergoing  a  rotation  and  translation  about  the  optical 
axis,  and  small  rotation  about  an  axis  orthogonal  to  the  optical  axis. 
More  generally,  if  h(x,t)  =  H(x,t)x,  where  H  is  a  3  x  3  matrix,  any 
rigid  body  motion  of  a  planar  section  in  the  scene  can  be  captured. 
Since  Q  can  be  approximated  by  planar  sections,  this  will  lead  to  a 
satisfactory  model  except  at  discontinuities  and  occluding  boundaries 
(relative  to  the  normal  to  the  image  plane).  Occlusions  are  represented 
by  a  real-valued  factor  m(x,t)  multiplying  the  right  hand  side  of  eq. 
28.  The  multiplier  m  may  depend  on  the  shape  (e.g.  curvature)  of  the 
surface,  illumination  of  the  surface,  background  lighting  and  possibly 
obstructions  in  the  held  of  view.  Clearly  0  <  m(x,t)  <  1,  where  1 
implies  that  the  point  x  on  the  surface  is  visible  at  time  f,  and  0  implies 
that  it  is  not.  Hence,  eq.  28  can  be  recast  in  the  form 

Ih,m(x,t)  =  m(x,t)I(h(x,t))  VT  G  D  x  [0,  oo)  fl/i-1(D)  C  M2x2,  (29) 

0  <  m(x,  t )  <  1. 

Since  Iihm{x ,  •)  must  evolve  continuously  in  time,  it  will  be  Marko¬ 
vian,  and  hence  St(I(x,t ))  :=  Ih,m{x,t)  =  m(x,t)I(h(x,t ))  is  a  one- 
parameter  semigroup  induced  by  the  how  h  :  Q  x  [0,  oo)  -G  M3  and  m  is 
a  cocycle  for  the  how  h,  i.e.  m  is  an  h-cocycle  Hence  the  image  deforma¬ 
tion  and  feature  tracking  within  a  given  view  can  be  stated  as:  Given  a 
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time  evolving  scene,  Ih,m  induced  by  the  flow  h  and  h-cocycle  m,  deter¬ 
mine  h  and  m.  This  problem  was  addressed  in  the  following  paper  [8]. 
In  this  paper,  the  class  of  all  semigroups  arising  from  cocycle-weighted 
time  warps  were  characterized.  Surprisingly,  it  was  shown  that  the 
cocycle-weighted  time  warps  give  rise  to  very  large  classes  of  flows  on 
Hardy  spaces  of  analytic  functions  and,  more  importantly,  a  relatively 
simple  criterion  was  developed  to  realize  the  weighted  deformation  as 
a  cocycle-weighted  time  warp.  This  simple  criterion  is  currently  be¬ 
ing  specialized  and  applied  to  answering  the  following  question:  Which 
filters  arise  from  time  warps?  [9]. 

Having  established  a  correspondence  between  the  multiple  camera 
views,  time  evolving  single  camera  views,  and  3-D  reconstruction  of  an 
RSO,  feature  detection  from  partial  data  is  of  utmost  importance  in  or¬ 
bit.  Here  the  basic  theory  for  processing  the  measurements  also  needs 
to  be  extended.  Standard  models  of  reconstruction  of  geometric  fea¬ 
tures  using  multiple  views  utilizes  data  in  each  view  and  backprojects 
this  data  to  construct  the  physical  object.  This  method  is  extremely 
costly  and,  in  the  long  run,  may  be  quite  inefficient  in  feature  detec¬ 
tion.  Thus,  it  was  asked:  Is  the  entire  data  in  each  image  necessary  for 
a  rapid  and  accurate  reconstruction  of  features  in  an  RSO?  This  ques¬ 
tion  led  us  to  investigate  spectral  methods  for  analyzing  the  camera 
images.  Description  of  the  2-D  views  in  appropriate  bases  that  remains 
invariant  under  the  image  formation  operator  (i.e.  the  /  operator)  al¬ 
lows  one  to  reduce  the  reconstruction  significantly,  while  extracting  the 
desired  features.  That  is,  we  construct  a  family  of  orthonormal  basis 
with  respect  an  appropriate  inner  product,  let  T  be  an  operator  that 
maps  the  2-D  image  into  £2  with  respect  to  this  basis  family,  and  re¬ 
quire  this  basis  to  have  the  property  that  T  commutes  with  the  image 
operator  I .  For  example,  if  T  is  the  Fourier  transform,  it  provides  a 
representation  in  terms  of  the  complex  exponential  basis,  and  T  com¬ 
mutes  with  a  slicing  operator.  Then  the  Fourier  transform  of  a  slice 
through  a  3-D  object  is  the  same  as  the  slice  through  the  3-D  Fourier 
transform  of  that  object.  This  is  the  projection-slice  theorem.  This 
theorem  allows  the  transformation  of  the  problem  from  a  pointwise  re- 


construction  of  a  geometric  object  from  multiple  views  to  a  spectral 
reconstruction.  Spectral  reconstruction  can  be  tailored  to  desired  res¬ 
olution  and  features.  Determination  of  how  the  transform  of  a  feature 
would  look  like  in  a  view  or  a  sequence  of  views  can  then  be  used, 
via  correlation  methods,  to  readily  detect  and  identify  those  features. 
The  results  of  this  partial  data  reconstruction  have  been  presented  at 
a  conference  [10]  and  are  being  prepared  for  publication  [11].  These 
results  establish  data  spaces  constructed  from  the  partial  data  which 
converge  in  weak  star  topology  to  the  complete  data  spaces.  Realiza¬ 
tion  of  specific  features  in  these  data  spaces,  allow  for  detection  of  the 
desired  features  at  much  smaller  cost.  A  closely  related  paper  to  this 
work,  dealing  with  representation  of  the  extreme  points  of  the  space 
of  positive  pluriharmonic  function  in  Euclidean  balls  was  also  devel¬ 
oped  [12]  under  the  AFOSR  grant.  The  geometric  tomography  method 
developed  this  paper  will  have  important  application  in  the  algorithm 
to  be  developed  to  describe  complex  images  from  multiple  views. 
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